Chargements des données et librairies

#installaation de package nécéssaire en local
#install.packages("BiocManager")
#BiocManager::install("phyloseq")
#install.packages("phyloseq")
#install.packages("remotes")
#remotes::install_github("mahendra-mariadassou/phyloseq-extended", ref = "dev")
#install.packages("gridExtra")


#Installation des librairies
library(tidyverse)
library(phyloseq)
library(phyloseq.extended)
library(ggplot2)    ## for graphs  necessary ???
library(ape)        ## for tree manipulation
library(vegan)      ## for community ecology analyses

library(scales)

library(grid)
library(gridExtra)  ## for design graph
library(knitr)
library(kableExtra)

chemin <- "~/save/Milint/Resultat/"

load("BIOM.Rdata")

Le BIOM a été créé dans le script01. Il est cependant brut

Etape de filtre sur l’abondance par échantillon et OTU

Somme de tous les reads et on applique un filtre de 10-5 afin de s’assurer que ce ne sont pas des artefacts Dans notre cas permet de filtrer les taxa avec des valeurs supérieures à 180 reads par OTU

filtre <- sum(taxa_sums(Biom))*1e-5

filter_biom <- prune_taxa(taxa_sums(Biom)>filtre,Biom)

save(Biom, filter_biom, file = "Data_Biom.Rdata")

Niveau d’abondance par échantillon

# abondance sur données complètes
sample_sums(Biom) %>% range()

[1]   5042 240312

# abondance sur données filtrées
sample_sums(filter_biom) %>% range()

[1]   5003 239715

# moyenne 
moy <- sum(taxa_sums(Biom))/846
moyf <- sum(taxa_sums(filter_biom))/846
mi <- min(taxa_sums(Biom))
ma <- max(taxa_sums(Biom))
mif <- min(taxa_sums(filter_biom))
maf <- max(taxa_sums(filter_biom))
med <- median(taxa_sums(Biom))
medf <- median(taxa_sums(filter_biom))

df <- data.frame(
  échantillon = c(846,846),
  mean = c(moy,moyf),
  median = c(med,medf),
  min = c(mi,mif),
  max = c(ma, maf)
)
rownames(df) <- c("Biom","Filter_Biom")
kable(df)

	échantillon	mean	median	min	max
Biom	846	21339.08	5	0	2709474
Filter_Biom	846	21207.16	1092	181	2709474

Le niveau d’abondance semble hétérogène. Il faudra donc prévoir un procédé de rarefaction

Graphique d’abondance des échantillons

données complètes

plot_bar(Biom)+
  ggtitle("Abondance par échantillon sur données totales")+
  theme(axis.text.x=element_blank(),
        axis.ticks.x=element_blank())

données filtrées

plot_bar(filter_biom)+
  ggtitle("Abondance par échantillon sur données filtrées")+
  theme(axis.text.x=element_blank(),
        axis.ticks.x=element_blank())

Dorénavant nous ne travaillerons que sur données filtrées Ne sert à rien d’aller au niveau Phylum ou infra car trop d’échantillons donc aucune couleur n’est vissible Nous remarquons que le niveau d’abondance par échantillon est hétérogène de 25 000 en moyenne à 250 000.

J’avais fait des niveaux d’abondance par échantillon en fonction de différentes variables (sexe, satut tabagique…) mais ne sert à rien car ça équivaut à voir la profondeur et ce n’est pas biologiquement informatif donc enlevé.

Répartition des OTU par phylum

t <- table(tax_table(filter_biom)[,"Phylum"])
kable(t)

Var1	Freq
Actinobacteria	22
Bacteroidetes	223
Firmicutes	686
Fusobacteria	2
Proteobacteria	45
Synergistetes	9
Tenericutes	5
Verrucomicrobia	5

Visualisation de l’abondance

Abondance en fonction des taxa

plot_bar(filter_biom, "Phylum", fill = "Phylum")

# Coloration ne fonctionne pas
 # scale_fill_manual( values = c("Actinobacteria" = "black", "Bacteroidetes" = "orange", "Firmicutes" = "blue",
 #                               Fusobacteria   = "yellow",
 #                                 Proteobacteria   = "red",
 #                                 Synergistetes= "pink",
 #                                 Tenericutes  = "green",
 #                                 Verrucomicrobia = "purple"))
 #  
  
#geom_bar(aes((x="Phylum) colour = "Phylum")

Bar plot abondance

#Abondance par phylum
phy <- tax_glom(filter_biom, taxrank = "Phylum")
plotdata<-psmelt(phy)
ggplot(plotdata,aes(y=Abundance))+
  geom_boxplot(aes(group="Phylum", fill="Phylum"))+
  facet_wrap(~Phylum, scales="free_y",nrow=2)

Abondance des firmicutes niveau famille au sein des échantillons

plot_composition(filter_biom, "Phylum", "Firmicutes", "Family", fill = "Family")+
  facet_grid(scales = "free_x") + 
  theme(axis.text.x = element_blank())

Abondance des firmicutes niveau des gènes au sein des échantillons Je ne vois pas l’interêt !!

plot_composition(filter_biom, "Phylum", "Firmicutes", "Genus", fill = "Genus", numberOfTaxa = 20) + 
  facet_grid(scales = "free_x", space = "free_x") + 
  theme(axis.text.x = element_blank())

Composition en Firmicutes au sein du groupe statut tabagique

plot_composition(filter_biom, "Phylum", "Firmicutes", "Genus", fill = "Genus", numberOfTaxa = 5) + 
  facet_grid(~tabac, scales = "free_x", space = "free_x")+
  theme(axis.text.x = element_blank())

Composition en Firmicutes au sein du groupe âge

plot_composition(filter_biom, "Phylum", "Firmicutes", "Genus", fill = "Genus", numberOfTaxa = 5) + 
  facet_grid(~AGE, scales = "free_x", space = "free_x")+
  theme(axis.text.x = element_blank())

* Composition en Firmicutes au sein du groupe sexe

plot_composition(filter_biom, "Phylum", "Firmicutes", "Genus", fill = "Genus", numberOfTaxa = 5) + 
  facet_grid(~SEX, scales = "free_x", space = "free_x")+
  theme(axis.text.x = element_blank())

Au vu des différents graphs, il y a une grande diversité et nous ne voyons pas de genre majoritaire quelque soit la variable. Pour l’ordre Clostrdiales est majoritaires quelque soit la variable prise (exemple ci dessous)

Composition en Firmicutes au sein du groupe âge

plot_composition(filter_biom, "Phylum", "Firmicutes","Order",  fill = "Order") + 
  facet_grid(~AGE, scales = "free_x", space = "free_x")+
  theme(axis.text.x = element_blank())

Je pense que les graphs suivants ne sont d’aucune utilité . J’ai fait les mêmes sur les femmes mais trop d’échantillons

knitr::opts_chunk$set(
  fig.width = 6, fig.height = 15, 
  fig.align = "center")

# échantillon ayant plus de 50000 abondance 
#prune_samples(sample_sums(biom) > 50000, biom)

#Abondance sur Hommes ayant plus de 25 000
H_more_25K <- subset_samples(filter_biom, SEX == 1 & sample_sums(filter_biom) > 25000)
H_m_25K <- plot_bar(H_more_25K, fill = "Phylum")

#Abondance sur Hommes ayant moins de 25 000
H_less_25K <- subset_samples(filter_biom, SEX == 1 & sample_sums(filter_biom) <= 25000)
H_l_25k <- plot_bar(H_less_25K, fill = "Phylum")
#plot_composition(H_less_20K, "Phylum", "Firmicutes", "Family", fill = "Family")

grid.arrange(H_m_25K,H_l_25k , nrow=2, 
             top=textGrob("abondance sup ou inf à 25K catégorie Hommes ",gp=gpar(fontsize=20,font=3)))

Alpha-diversité

filter_rare <- rarefy_even_depth(filter_biom, rngseed = 20210416)

Graphique Covariable et test statistique

Exploration de l’impact de chaque covariable sur la diversité-alpha

plot_richness(filter_rare, x = "AGE", measures = c("Observed", "Shannon", "InvSimpson")) + 
  geom_boxplot(aes(group = AGE)) ## add one boxplot per AGE  values

# Couleur ne fonctionne pas scale_fill_manual(values=c("orange", "blue", "green3","red","yellow"))

Age diffère significativement en terme d’observation d’OTU

div_data <- cbind(estimate_richness(filter_biom),  ## diversity indices
                  sample_data(filter_biom)         ## covariates
                  )

model <- aov(Observed ~ AGE, data = div_data)
anova(model)

Mais également en terme de diversité en nombre d’espèce

model <- aov(InvSimpson ~ AGE, data = div_data)
anova(model)

model <- aov(Shannon ~ AGE, data = div_data)
anova(model)

plot_richness(filter_rare, x = "SEX", measures = c("Observed", "Shannon", "InvSimpson")) + 
  geom_boxplot(aes(group = SEX)) ## add one boxplot per SEX  values

Pas de différence signicative pour le nbre d’OTU mais différence significative pour le nbre d’espèces

model <- aov(Observed ~ SEX, data = div_data)
anova(model)

model <- aov(InvSimpson ~ SEX, data = div_data)
anova(model)

#model <- aov(Shannon ~ SEX, data = div_data)
#anova(model)

plot_richness(filter_rare, x = "tabac", measures = c("Observed", "Shannon", "InvSimpson")) + 
  geom_boxplot(aes(group = tabac)) ## add one boxplot per AGE  values

Analyse Biome

DUBii 2021

16 avril, 2021